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Abstract 


Pulsars are a form of variable celestial source that have shown to be usable as aids for autonomous, deep space 
navigation. Particulary those sources emitting in the X-ray band are ideal for navigation due to smaller detector 
sizes. In this paper X-ray photons arriving from a pulsar are modeled as a non-homogeneous Poisson process. 
The method of pulse phase tracking is then investigated as a technique to measure the radial distance traveled 
by a spacecraft over an observation interval. A maximum-likelihood phase estimator (MLE) is used for the case 
where the observed frequency signal is constant. For the varying signal frequency case, an algorithm is used in 
which the observation window is broken up into smaller blocks over which an MLE is used. The outputs of this 
phase estimation process were then looped through a digital phase-locked loop (DPLL) in order to reduce the 
errors and produce estimates of the doppler frequency. These phase tracking algorithms were tested both in a 
computer simulation environment and using the NASA Goddard Spaceflight Center X-ray Navigation Laboratory 
Testbed (GXLT). This provided an experimental validation with photons being emitted by a modulated X-ray 
source and detected by a silicon-drift detector. Models of the Crab pulsar and the pulsar B 182 1-24 were used 
in order to generate test scenarios. Three different simulated detector trajectories were used to be tracked by 
the phase tracking algorithm: a stationary case, one with constant velocity, and one with constant acceleration. 
All three were performed in one-dimension along the line of sight to the pulsar. The first two had a constant 
signal frequency and the third had a time varying frequency. All of the constant frequency cases were processed 
using the MLE, and it was shown that they tracked the initial phase within 0.15% for the simulations and 2.5% 
in the experiments, based on an average of ten runs. The MLE-DPLL cascade version of the phase tracking 
algorithm was used in the varying frequency case. This resulted in tracking of the phase and frequency by the 
DPLL outputs in both the simulation and experimental environments. The crab pulsar was experimentally tested 
with a trajectory with a higher acceleration. In this case the phase error tended toward zero as the observation 
extended to 250 seconds and the doppler frequency error tended to zero in under 100 seconds. 
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1 Introduction 


Throughout human history, celestial sources have been used as aids for navigation. Most of these sources are 
constant in emissions, and serve as optical point source references, whose use dates back to ancient sailers navigating 
by the stars. A small portion of celestial sources emit radiation at varying intensity. This variation allows for new 
possibilities for use in navigation algorithms that were not acheivable with steady point sources. There are many 
causes of these varying emissions including object pulsations, eruptions, rotation, eclipses by companion stars, and 
cataclysmic effects [1] [2]. A subset of these varying celestial sources are spinning neutron stars or pulsars. They 
were first discovered in the radio band in 1967 [3]. Subsequently, pulsars have been discovered that emit radiation 
throughout the electromagnetic spectrum [4] [5] . When it comes to navigation, those that emit in the X-ray band 
(0.1 - 100 keV) are the most attractive because smaller detectors can be utilized than for optical or radio sources. 
Figure 1 is an image of the Crab Pulsar in the optical and X-ray bands. 

Neutron stars are formed as the result of a supernova explosion of a massive star that has run out of fuel [6] 
[7]. Conservation of angular momentum spins these stars up to high rates during their collapse. Neutron stars 
also have very strong magnetic fields, which accelerates charged particles along their field lines [8]. This results 
in powerful beams of electromagnetic waves, across the entire spectrum, being radiated from the pulsar’s magnetic 
poles. Misalignment between the spin axis and magnetic field axis causes an observer to be able to sense a pulse of 
magnetic radiation due to the magnetic pole sweeping across the line of sight [8]. 

Pulsars have very regular, stable and periodic signals [4] [5]. In fact they are so stable that some have been shown 
to match the stabilities of today’s atomic clocks [9] [10]. Pulsars have unique periodic pulse profiles, which along with 
their relatively uniform distribution throughout the sky, suggests they might be able to be used as a natural celestial 
lighthouse system for navigaion. The baseline X-ray pulsar-based navigation (XNAV) method uses observations of 
millisecond period pulsars (MSPs). 

Due to their great distance and dispersion, pulsars and XNAV specifically are a method for autonomous navigation 
in deep space. Some current methods of spacecraft navigation throughout the solar system include: vehicle tracking 
using the Deep Space Network (DSN) [11] , occultation of celestial objects [12], and orbit-determination using 
Earth-based observations. The primary method is DSN radiometric tracking. This provides a range and range-rate 
measument along the line-of-sight from the Earth to the spacecraft. Differential one-way ranging (Delta- DOR) 
using very long baseline interferometry (VLBI) with the DSN can acheive angular measurements with accuracies 
on the order of 1 nrad [13]. XNAV measurements have the same accuracy throughout the solar system, whereas 
measurements from the DSN decrease in accuracy with increasing distance from Earth. Thus, for distant missions, 
XNAV measurements can either be added to DSN measurements to increase accuracy, or be used alone in order to 
eliminate the reliance on an uplink from Earth. 

Most of the astronomical observation of pulsar sources have been done in the radio band because these emissions 
can be received on the Earth’s surface. This has resulted in suggested methods to use pulsar radio signals for 
navigatioin [14], but significant improvement would need to be made in detector technology in order to make this 
feasible for space missions. Consequently, there is greater focus on methods to observe pulsars in the X-ray band to 
test a possible XNAV system and provide updated timing models. The All-Sky Monitor (ASM) instrument aboard 
the Rossi Timing Explorer (RXTE) has served the role of monitorying the X-ray sky in the past [15]. In the future 
the Neutron-star Interior Composition ExploreR (NICER) and Station Explorer for X-ray Timing and Navigation 
Technology (SEXTANT), in development at the Goddard Space Flight Center (GSFC), will serve these purposes 
attached to the International Space Station (ISS) [16] [17]. These combined missions will perform phase-resolved 
spectroscopy of rapidly spinning neutron stars and demonstrate real-time X-ray pulsar navigation using MSPs, 
respectively. 

1.1 Background 

A pulsar’s timing model can be accurately measured over long observations. These models can be represented as 
the total accumulated phase of a pulsar’s signal over time. An initial signal number <l>o is chosen at a fiducial time, 
to at the Solar System Barycenter (SSB), which is used as an inertial reference location. The total phase can be 
modeled as a function of time and broken down into fractional and whole integer components [8] , 

m = w) + N(t), (i) 

where <f> is the fractional portion of the pulse period and N is the number of integer cycles since to- This can also be 
expressed in terms of the pulse signal frequency and its derivatives. This is know as the pulsar spin equation [8], 

$(t) = $(t 0 ) + f[t — to] + ^[t — to] 2 + ^[t — to] 3 - (2) 
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Figure 1: On the Left is an image of the Crab Nebula and Pulsar, in the optical (red) and X-ray (blue) portions 
of the spectrum superimposed. The right image is the Crab Pulsar in only the X-ray potion of the spectrum 
(NASA / CXC / ASU / J . Hester et al.) [18], 


In order to be used in navigation, pulsars need to be accurately modeled and identified. Previous methods of 
determining a pulsar source’s characteristics have relyed on creating binned pulse profiles over long observations 
of a source and then folding back over the expected pulse period [1]. Observed pulse time of arrivals (TOAs) are 
then computed through comparison with a high signal-to-noise ratio (SNR) template. The pulse signal needs to be 
accurately timed over each observation with respect to an inertial system. This time transformation must be done 
carefully to avoid a smearing effect that would be present if removal of the motion and higher order Doppler effects 
is done incorrectly, causing increased navigation errors [1]. Various methods for postion determination have been 
studied in the literature [19], [20], [21], [22] [23], and [24], These approaches fall in two general categories: absolute 
3-dimensional position and velocity determination in an inertial frame, and delta-correction where estimates of range 
and range-rate values are updated. Both methods provide the means to maintain a continuous navigation solution 
[ 8 ], 


1.2 Suggested Approach 

The approach used in this paper is modeled after the pulse phase tracking algorithm that was developed in [1]. 
This method attempts to exploit the periodic nature of pulsars by measuring the pulse TOAs onboard a moving 
spacecraft. The pulse phase measured on the spacecraft depends on the radial distance between the spacecraft 
and the pulsar. Thus, information about the spacecraft’s position can be extracted from these measurements. The 
practicality of using this for navigation purposes was looked at in the radio band in [19] [25] [26] and in the X-ray 
band in [24] [25] [27], 

The method of phase tracking relies on the use of a maximum- likelihood phase estimator (MLE), which can 
estimate the initial phase at the beginning of an observation of a pulsar given a set of TOAs over that interval. 
When the signal frequency is constant at the spacecraft this measurement of the initial phase is enough to know the 
phase evolution over the entire observation time interval since it changes linearly with known slope. In the case of 
changing signal frequency, due to the dynamical motion of the spacecraft, the observation window can be broken up 
into smaller blocks over which the frequency can be assumed constant. Then the phase of the observed pulse can be 
tracked over the observation as the spacecraft travels along its trajecoty, which gives a measure of distance and range 
rate [1]. The output of the MLE in this case is sent through a digital phase-locked loop filter (DPLL) at each iteration 
in order to reduce the errors. An advantage of this method is that, if it is used as a delta-correction technique, it can 
possibly provide continuously updating position information as opposed to the infrequent TOA-difference technique 
[ 8 ]- 

The proposed method tracks the expected parameters of the source signal in order to estimate the spacecraft 
trajectory. This allows continuous updates of the vehicle motion over short time frames. An added advantage of 
phase tracking is that it eliminates the requirement of time transferring all of the photons arriving at the spacecraft 
to an inertial reference at the SSB before measurements can takes place. This eliminates the risk of smearing the 
phase profile over the time transfer. This phase tracking algorithm was demonstrated through simulation in [1] by 
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looking at an RXTE observation of the Crab pulsar. The algorithm was shown to have good phase and frequency 
tracking capabilities. The goal of this paper is to experimentally validate this phase tracking algorithm using a 
lab setup that includes a modulated X-ray source to emit photons that model those arriving from a pulsar and a 
detector to read and time-tag the collection of the emited X-ray photons. This setup was then used to perform phase 
tracking on the received TOAs at the detector in the lab. This adds the complexity of dealing with actual photons 
and detector physics outside of the simulation environment. This method was tested over various scenarios over a 
range of constant and changing signal frequency cases. Both the Crab pulsar (B0531+21) and B1821-24 were used 
as source models and the results were compared to simulations. 

2 X-Ray Pulsar Signal Model 

2.1 Non-Homogeneous Poisson Process 

Due to the periodic nature of pulsar signals, the time of arrivals of X-ray photons at a spacecraft detector can be 
modeled as a non-homogeneous Poisson process (NHPP)[1]. With this assumption, the number of photons arriving 
from a pulsar in a given interval is a Poisson random variable with a periodic or possibly quasi-periodic rate function 
A (t) > 0. Given this model, the probability that k photons arrive from a source in an interval (a,b) is given by the 
following [1] [28], 

Pr[fc; (a, 6)] = ^y{exp[— J \{t)dt]}[j A {t)dt] k . (3) 

The rate function A (t) is a cumulation of all photons that hit the detector. This includes both the source photons, 
which arrive at an average rate of R s photons per second (ph/s) and the background photons arriving at the constant 
rate R b ph/s [1]. The source photons include both pulsed and steady emissions from the pulsar, and the pulse fraction 
p gives their ratio. The non-pulsed source photons behave like background photons, which makes it difficult to tell 
them apart at the spacecraft detector. The following consequences stem from the use of the NHPP model [1]: 

1. Setting k = 0 in Equation (3) gives the probability that zero photons will arrive in (a,b) 

Pr[0; (a, b )] = exp[— f A {t)dt] (4) 

J a 

2. The Probability that exaclty one photon arrives in an infinitesimal interval of length At centered at time t is 
given by evaluating Equation (3) with k = 1 

Pr[l; (-At/2, +At/2)] = A(t;6» 0 ,/)At as At 0, (5) 

where A(t;0o;/) is the rate function given in terms of the initial phase and frequency parameters, which will be 
discussed later. 

3. The Probability of more than one photon arriving in an infinitesimal interval of length At is zero when At — > 0. 

4. The number of photons arriving in a given interval is independent from those of all other disjoint intervals. 

2.2 Pulsar Rate Function 

The time- varying rate of arrival of the pulsed source photons can be modeled as the composition of the normalized 
pulse profile function, h{(/>), and observed phase at the spacecraft detector, (j>det{t ) [!]■ This, along with the notion 
of constant rate background and non-pulsed source photons, allows us to write the overall rate function for a given 
pulsar as [1] 

\{t) = R b + (1 - p)R s + pR s h(4>det{t)) {ph/s). (6) 

We can define parameters a and /3, called the effective source and background photon arrival rates, as follows 


a = pR s {ph/s) 

P = R b + (1 - p)R s {ph/s). 


(7) 
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These parameters capture the flux of the periodic source photons and the constant flux of the cumulation of the 
background and non-pulsed source photons, respcectively. They are important for considering performance, and 
their ratio corresponds to the signal-to-noise ratio (SNR) of the pulsar with a particular detector setup. 

The pulse profile function h(<j)) is normalized to satisfy the conditions: 


1 . min^ 0tl )h((j>) = 0 

2 . £ h(4>)d(j) = 1 . 

Usually, h is defined for <j> G (0,1), but in order to be used in Equation (6) it needs to be extended as a peri- 
odic function with a period of one cycle. This can be done by letting h(m + </>) = Vto e Z [1], Then h is 
defined for <j> G (— oo,+oo), which allows h(<j>det{t )) to continue to be computed in Equation (6) even as the phase 
accumulated at the detector grows over the observation interval. 

Now, define the phase observed at the detector as the sum of the initial phase, do, and the accumulation of phase 
since the start of the observation period, 6(t). Then in terms of the observed signal frequency, /(£), the accumulated 
phase at the detector over an observation interval (to,t) is given as 


<t>det{t) = 0o + 9(t) — do + f (9) 

Jt 0 

The observed signal frequency can be broken up into source frequency and doppler frequency componenets [1]: 

f(t) = fs+ fd(t) = f s ( 1 + v(t)/c), (10) 

where the doppler frequency shift is fd = f s v(t)/c , v(t) is the range rate of the spacecraft in the line-of-sight direction 
to the pulsar, and c is the speed of light. 

For the sake of this paper, the source frequency, f s , is taken as constant on the observation interval, and the 
second and higher-order relativistic doppler effects are ignored. Under these assumptions, the pulsar spin equation, 
Equation (2), simplies to Equation (9). This will allow for the development of a Maximum-Likelihood Estimator and 
Phase- Tracking Algorithm as in [1] . As explained in [1] , the scope of analysis is not limited by these assumptions since 
source frequency variations are often forcasted, which allows for pre-processing of the photon TOA measurements to 
remove the effects due to these changes. 

Under these assumptions, the phase accumulated since the beginning of the observation has two terms [1], 

0(t) = f a {t - t 0 ) + [ fd(t')dt', (11) 

Jt 0 

where 0d{t) = f* fd(t')dt' is the doppler phase. If range-rate v(t) is a constant v, 0(t) = / s (l + v/c)(t — t 0 ) is linear 
and thus the Poisson rate function A (t) is periodic. For non-constant range-rate, the Poisson rate function is only 
quasi-periodic due to the nonlinear dependence of 0d(t) on time. The cases of constant and time-varying observed 
signal frequency can be considered and dealt with separately. Now look at the rate function under each scenario. 

2.2.1 Constant Observed Signal Frequency 

First, notice that with constant observed signal frequency, / = f s ( 1 + v/c), the phase at the detector is given 
linearly as follows [1], 

<Pdet = $o + f(t — to)- (12) 

Solve for the Poisson rate function by substituting Equation (12) into Equation (6), which gives 

■Mt; 00) /) = P + oth{do + f(t — to))- (13) 

This shows the functional dependence of the rate function on do and / and that A is strictly-periodic due to the 
periodicity of h. When the detector is stationary or is moving at a constant range-rate with respect to the pulsar, 
the constant signal frequency model can be used. Also, the period change of the pulsar over the observation time 
interval should be negligible[l]. 
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2.2.2 Time- Varying Observed Signal Frequency 


The observed phase at the detector has a term contributed by the doppler effects when the frequency is varying 
[1], In this case, 

<t>det = $0 + f = 9 q + f [f s + fd(t')\dt' = 0q + f s (t — to) + 9 d(t) . ( 14 ) 

Jt 0 Jto 

Now plug this into Equation (6) to get the quasi-periodic Poisson rate function for the time-varying frequency case 

X(t) = p + ah(6o + f s (t — to) + 9d{t)). ( 15 ) 

This model must be used if the radial speed with respect to the pulsar is varying considerably over the observation 
interval. 


3 Maximum-Likelihood Initial Phase and Frequency Estimation 

Now the Maximum-Likelihood Estimator (MLE) for 9 0 and f can be derived following the approach of [29] and 
[1] for the constant frequency model. The goal is to estimate these parameters based on a set of TOAs, which are 
either a computer simulated realization of the NHPP or an experimentally detected realization of the NHPP. For 
background on MLEs [29], [30], and[31] were consulted. 

Assume that the pulse profile function h(<j>), effective source arrival rate a, and effective background arrival 
rate /3 are all known. Consider an observation interval (t 0 ,t 0 + T obs ) and a set of photon TOAs {t k }// =0 with 
t 0 <t k < to + T obs Vfc. These TOAs are a realization of the NHPP for the pulsar with rate function A (t; 6 0 , f) given 
by Equation (13). The aim is to be estimate the model parameters 9q and /. 

The first step is to compute the K-dimensional joint probability density function (pdf) of the NHPP realization, 
/({tfc}). This is acheived by considering infinitesimal intervals of length A t k centered around each t k , k = 1, 2, • • • , K. 
The joint pdf can be computed by first looking at the probability that one and only one photon arrives within each 
interval and none outside of them [1] 

f({t k }) At! At 2 • • • A t K = Pr[0; {t 0 ,h - Ah/2)} x Pr[l; (h - Ah/2, h + Ah/2)} (16) 

xPr[0; (h+ Ah/2, h-Ah/2))} x ••• x Pr[l; (t K - At*/2, t K + At K /2)} 
xPr[0; {t-K + At K /2,to +T obs )}. 

From the consequences of the NHPP model plug in Equations (4) and (5). Now take A t k 0 for all 1 < k < K. 
This yields, 


r~to-\-T 0 bs 

f({tk}) = exp [ - / A (t; 9 0 , f)dt ] A (t k ; 9 0 , f), (17) 

Jt ° k = 1 

which is the likelihood function [1], The desired MLE solves for the values of 9 0 and / at which the likelihood 
function is maximized. It is equivalent to consider the logarithm of the likelihood function (LLF), which is easier to 
minimize, 

~ ~ ^ „ rtod-Tobs ~ _ 

LLF(0 0 , f) = ^2 log[A(< fc ; 9 0 , /)] - / A (t k \ 0 O , f)dt. (18) 

fc= l ^ to 

Here log is the natural logarithm, and do and / are the values at which the LLF is computed. In most cases, the 
integral in Equation (18) can be dropped since it displays a minimal dependence on the search model parameters. This 
minimal dependence holds when the observation interval lasts many cycles of the signal at the expected frequency 
( fT 0 bs ^ 1) [1]- The optimization problem that needs to be solved by the MLE is: 

K 

0o, /) = argmax Soee j e n ^ log[/3 + ah{9 0 + f(t k - t 0 ))], (19) 

k — 1 

where 0 and fl are the search spaces for phase and frequency, respectively. Solving this maximization by brute-force 
involves nested, iterative grid-searches. Each iteration searches a smaller region, but with greater precision around 
the output of the previous iteration. The simulation methods used specifically for this paper will be discussed in a 
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later section. There are also more sophisticated methods that attempt to decrease the MLE run time. For example 
[32] uses a Discrete Fourier Transform procedure to increase speed, but it sacrifices the gaurantee that the output is 
a maximum. 

If it is assumed that the signal frequency and doppler shift are known, then Equation (19) reduces to a phase-only 
optimization, 

9 0 = argmaXg oee LLF(e 0 , /„/), (20) 

where / re / is the a-priori known signal frequency. It was shown in [1] that the performance of this MLE is optimal 
since it acheives the Cramer-Rao performance bound after a certain threshold. If the observed signal frequency is 
constant, the measurement of the initial phase 9q by the MLE and knowledge of the frequency /, either a-priori or 
from the MLE, is sufficient to determine the phase over the entire observation interval. This is because the phase at 
the detector has a linear dependence on 9q and / if the signal frequency is constant, as seen in Equation (12). 

4 Pulse Phase Tracking 

When observations are made onboard a spacecraft whose trajectory curves with respcect to the pulsar, or the 
radial speed in the direction of the source varies by a noticable ammount, then the observed frequency can no longer 
be considered constant without acruing errors that will overwhelm the navigation scheme. The procedure of pulse 
phase tracking of [1] will be adopted to deal with this. This is based on the time-varying signal frequency model 
introduced in Section 1.2.2. The first step is to partition the observation interval into N blocks of size T, which is 
set small enough so that the observed signal frequency is approximately constant over each block. This is illustrated 
in Figure 2. Photons arriving in the n-th block are processed in a batch using the MLE of Equation (20) to solve for 
the y-intercept 9o(n) as seen in Figure 2. When the reference frequency is equal to the source frequency, 9q(ji) is an 
estimate of 9q + 9d(t) for t = nT. Thus, by using MLEs over one block of photon TOAs at a time, this algorithm 
tracks the phase quantity 6*o + 9d(t) over the entire obseration interval [1]. 

There are a few issues that need to be addressed with regards to implementing this phase tracking algorithm. 
First, the sequence of estimates 9o(n) may contain jumps due to phase wraparounds. In order to deal with this, 
the search space for the first block will be the entire phase space 0 = (0,1), but subsequent blocks will take the 
previous block’s phase estimate 9 0 {n — 1) as the center of the search window. This allows the sequence 9o(n) to 
wander outside of (0,1), if necessary, to avoid jumps. The other main issue is the presence of noise due to unmodeled 
vehicle dynamics and estimation errors. To reduce these effects, the sequence of outputs 9o(n) from the MLE can be 
post-processed through a digital phase-locked loop (DPLL). This allows for a sequence of estimates 9 EPLL {n) and 
fd PLL ( n ) ^0 track the phase and doppler frequency with smaller errors than using the MLE alone. 

4.1 Digital Phase-Locked Loops 

A digital phase-locked loop (DPLL) filter can be used to reduce the effects of the inherent estimation errors and 
noise in the MLE phase outputs #o(n) [1]. Phase-Locked Loops are abundant in modern signal processessing and 
communication systems [33]. The idea is to be able to form a local representation that is synchronized with an 
incoming signal’s frequency and phase. An overview of phase- locked loop techniques is included in [34]. Originally, 
all phase-locked loops were analog, and their implimentation and design relyed extensively on knowlegde of analog 
systems [35]. Techniques for designing and implementing analog phase-locked loops are discussed in [36] and [37]. 
In order to take advantage of advances in digital computing, fully digital loops have been developed. DPLLs have 
increased in popularity due to andvantages with respect to performance, speed, reliability, size, and cost [33]. Also, 
DPLLs allow for their parameters to be defined such that each has a direct physical meaning with respect to loop 
noise bandwidth Bl, root-specific decay rate, or root-specific damping [38]. Digital Phase-Locked loop design and 
implimentation are discussed at length in [33]. 

The second-order DPLL used in [1] and developed in [38] will be presented and used in the pulse phase tracking 
algorithm. This formulation allows the loop-filter constants to be determined from loop roots that can be placed 
anywhere in the complex plane [38]. The DPLL will be implemented in a cascade with an MLE as seen in Figure 3. 
Thus, it is assumed that the photon TOAs arriving at the detector are first processed through the MLE, which gives 
an estimate 9Q ILE (n) at each block. The reference frequency used in the MLE is set equal to the source frequency, 
which is assumed constant, f re f = f s . Now the algorithm for the second-order DPLL will be presented. The first 
step is to compute the difference between the input phase from the MLE and the DPLL’s phase estimate. This is 
done through the following ideal phase detector, 

4>(ri) = 0§ ILE (n) — 9q PLL (ti). ( 21 ) 
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Figure 2: This diagram demonstrates the pulse phase tracking concept presented in [1], The observation interval is 
broken up into N blocks of size T seconds over which the signal frequency is approximately constant. 


This is carried out without wrapping around the phase. Next, the state of the loop filter can be updated using the 
error feedback (j>(n ) [38] 

n 

fd PLL (n + 1 )T = Ki4>(n) + if 2 ]T (22) 

m= 1 

where f PPLL is the measurement of the doppler phase by the DPLL and T is the loop update interval, which is set 
equal to the block duration of the phase tracking algorithm. This sequence, f EPLL {n + 1 )T, drives a numerically 
controlled oscillator (NCO), which sets the amount the current NCO phase must be changed in order to achieve the 
phase estimate of the next block. The NCO corresponds to the following phase accumulation at baseband, 

0o PLL (n + 1) = 0 PPLL {n ) + f PPLL (n + l)T. (23) 

In order to initialize the phase and frequency states of the DPLL before the start of the loop take 0 PPLL ( 0) = 9q ILE (0) 
and f PPLL ( 0) = f 3 v(to)/c. If there is not an available estimate of the initial range rate v(t 0 ) set f PPLL ( 0) = 0. 

5 Simulation Structure 

First, the MLE and pulse phase tracking MLE-DPLL cascade were tested in a computer simulation environment 
using MATLAB. In this case the photon TOAs received by a detector are generated by a realization of the NHPP 
with rate function given by Equation (6). This set of TOAs is first generated as if they are located at the solar 
system barycenter (SSB) and then delay is added in order to acount for the offset between the detector and the SSB. 
The NHPP photon generation is completed using the inversion method of non-uniform random variate generation 
presented in [39] and implemented for the specific example of X-ray pulsar photon generation in [40] . The general 
relativistic delay and doppler effects must be taken into account to tranfer the barycentric set of time of arrivals to 
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Figure 3: This figure shows the MLE-DPLL cascade used to reduce errors in the pulse phase tracking algorithm 
when the observed signal frequency is time- varying [1]. 


a set arriving at the spacecraft detector, which is desired to test the algorithms. To first order, this time transfer is 
given by 

tsSB - t sc = , (24) 

c 

where tssB is the time at the SSB, t sc is the time at the spacecraft, n is the unit normal vector pointing from the 
SSB to the pulsar, and r is the positon vector of the spacecraft with respect to the SSB, as seen in Figure 4 [41]. 
This set of TOAs at the spacecraft detector is what is sent to the MLE or MLE-DPLL cascade algorithms. 

In the constant signal frequency case the MLE algortithm is run one time using all of the TOAs in the entire 
observation interval. This results in estimates for the phase at the detector at all times in the observation, because 
it is given linearly by Equation(12) and the MLE gives estimates of the unknown model parameters. This method 
is used to test both cases where the detector is stationary and where the detector is moving with a constant range 
rate with respect to the pulsar. When the signal frequency is time-varying the observation interval is broken up into 
blocks and the pulse phase tracking algorithm is used along with the MLE-DPLL cascade structure. 



Figure 4: This figure shows the relationship between the spacecraft position, the SSB, and the pulses arriving from 
the distant pulsar, which can be assumed to be arriving as uniform wave fronts throughout the solar system [41]. 


6 Experimental Procedure 

The main goal of this paper is to experimentally validate the phase-tracking algorithm. The lab equipment of 
the NICER and SEXTANT teams at the NASA Goddard Space Flight Center (GSFC) was used for this purpose. 
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The GSFC X-ray Navigation Laboratory Testbed (GXLT), that the SEXTANT mission team used for their own 
development, was ideal for testing the algorithms discussed in this paper. The GXLT is made up of four major 
components: the X-Ray Digital Pulse Generator (XDPG), which creates a real-time implementation of the pulsar 
rate signal given in Equation (6), the Modulated X-Ray Source (MXS), which emitts the X-ray photons, the Silicon 
Drift Detector (SDD), which receives the photons and produces a single short pulse using a pulse shaping filter, 
and the Time- Tagging X-Ray Detector (XDET), which puts a time stamp on the photon pulses from the SDD [17]. 
The XDPG takes the pulsar definiton and the desired phase dynamics and creates a real-time implimentation of 
the rate-function at the spacecraft and then drives the input of the MXS as an NHPP [17]. Figure 5 illustrates the 
interconnections of the lab equipment and their use in generating the desired X-ray photons. 

The MXS is discussed in [42], and it uses an Ultraviolet (UV) Light Emitting Diode (LED) to modulate, through 
photoelectron flux, the X-ray photon emissions [17]. This choice of source allows for better modeling of MSP X-ray 
photons because it allows for fast switching response, intensity control, and is completely off when switched off [42]. 
Other common X-ray sources lack some of these requirements. Analysis of the SDD is included in [43]. The choice 
to use an SDD by the NICER team was due to it pocessing excellent energy resolution, operational simplicity, and 
fairly good time resolution [43]. Figure 6 contains a picture of the lab setup during a run. The oscilloscope shows 
the shape of the desired pulse profile as well as an output of the photons as they arrive at the detector. The MXS 
is commanded to emit photons that model TOAs received at the desired spacecraft location. This produces a set 
of TOAs received by the SDD that corresponds to photons the detector would see if it was traveling on the desired 
trajectory. The output of the SDD is sent through a pulse shaping filter, which produces a Gaussian pulse-like 
voltage signal [17]. This signal can be observed on an oscilloscope, as in Figure 7, in order to see the arrival rate 
matches the XDPG driving signal. 

Although this procedure uses artificially generated photons from an MXS instead of actual detection of a pulsar 
signal in the X-ray band, it still provides the added complexity and noise related to generating actual X-rays and 
detecting them with an SDD. This is a step above the purely computer simulated results because it includes the 
processing of X-ray TOAs at a detector, which was unaccounted for in the simulation. The output no longer 
automatically matches the realization of the NHPP. The same test situations for the computer simulations were used 
with the experimental equipment, and the results were compared. 


Vacuum Chamber 




Figure 5: This is a diagram demonstrating how the XDPG, MXS, SDD, and XDET work together to generate and 
detect the desired X-ray photons. The internal workings of the X-ray generation process performed by the MXS is 
included. The oscilloscope output image is taken from [17]. 
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Figure 6: This is a picture of the lab setup for the X-ray navigation ground testbed in use by the NICER team at 
GSFC. Image taken from [17]. 



Figure 7: This image shows the XDPG driving signal in red and the corresponding photon arrival events at the SDD 
for a low SNR. Image taken from [17]. 


7 Results and Discussion 

7.1 Simulation Results 

The pulsars that were used to test the phase tracking algorithms with simulations were B0531+21, the Crab 
pulsar, and B1821-24. The parameters for each of these pulsars that were used for the simulations are given in Table 
1. An amplified version of B1821-24 needed to be used to implement the phase tracking by blocks algorithm, since 
otherwise the flux was too low for a sufficient number of photons to arrive in each block. This artificial amplification 
is for the purpose of demonstrating the validity of the algorithm and attempts at a more realistic implementation, 
over a wider range of pulsars, will be investigated in future work. The amplified version of B1821-24 was also used 
for the constant frequency MLE cases. This was in order to make the observation times be equal to those used for 
the Crab pulsar. Figures 8 and 9 contain the pulsar pulse profile templates that were used to model the Crab pulsar 
and B1821-24, respectively. 
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Table 1: Pulsar Parameters for Simulation Runs 


Pulsar fs (Hz) a (ph/s) 

P (ph/s) 

B0531+21 (Crab) 
B1821-24 
B 182 1-24 amplified 

29.9401 

327.86885 

327.86885 

1078 

0.18914 

1470 

462 

0.00386 

30 


PSR B0531+21 (Crab) Data File 



Figure 8: This is the template for the Crab pulsar profile used in the simulations, shown over two periods. 


PSR B1821-24 Data File 



Figure 9: This is the template for the pulse profie of B1821-24 used in the simulations, shown over two periods. 


Simulations were developed around three basic scenarios in order to test the phase tracking algorithms. Two fall 
under the constant signal frequency model and the third requires the signal frequency to be considered time- varying. 
All three assume a one-dinrensional orbit determination problem, the spacecraft is moving along the line-of-sight 
direction to the pulsar. This simplificaiton allows the orbit information for the spacecraft to be easily determined 
and reduces the complexity of the problem, so it is easier to focus on the performance of the phase tracking algorithms. 
Future work will expand this to three-dimensional orbits. The two constant signal frequency scenarios considered 
were a stationary detector and a detector moving at a constant velocity. For the time- varying frequency scenario, the 
spacecraft is taken to be moving with a constant acceleration. These cases all allow for easy analytical determination 
of the truth trajectories for the phase at the detector. For instance, we know that in the constant signal frequency 
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case, the phase evolves linearly at the detector and thus the truth trajectory is determined completely by the initial 
phase 0o and the observed signal frequency / = f s (l + v/c). For the time- varying signal frequency case with constant 
acceleration a, the doppler frequency shift is given by 


fdit) = f s 


Vq + at 
c 


This can then be integrated in Equation (14) to give the phase at the detector 


(25) 


<f>det(t) =9 0 + f s (t - to) + ^[v 0 (t - to) + |(t - to) 2 ]- (26) 

The phase tracking algorithm provides estimates of (f>det(t) — fs{t — to), so the truth trajectory for the constant 
acceleration case is 

0o(t) = 9 0 + ^[v 0 (t - to) + |(t - to) 2 ]- (27) 

Table 2 lists all of the different simulations that were run along with their parameters. For the constant signal 
frequency cases the TOAs were processed through the MLE of Equation (20) in order to give estimates 9o ■ For each 
simulation Monte-Carlo methods were used on ten repeated trials in order to attempt to remove some of the affects 
due to non-deterministic effects, such as noise, on the accuracy of the MLE algorithm. The results for the stationary 
and constant velocity cases are presented in Table 3. The simulated results all had less than 0.15 % error from the 
true initial phase of 9 0 = 0.5. An example of an MLE plot for one of the scenarios is provided in Figure 10. 


Table 2: Descriptions of the Performed Simulation Runs 


Pulsar 

Type of Run 

00 

Vo (km/s) 

a (km/s 2 ) 

Initial fd (hz) 

Number of Runs 

Observation Time (sec) 

Crab 

Stationary 

0.5 

0 

0 

29.9401 

10 

10 

B1821-24 

Stationary 

0.5 

0 

0 

327.86885 

10 

10 

Crab 

Constant v 

0.5 

0.1 

0 

29.9402 

10 

10 

B1821-24 

Constant v 

0.5 

0.1 

0 

327.86995 

10 

10 

Crab 

Constant a 

0.3 

0 

0.1 

29.9401 

1 

200 

B1821-24 

Constant a 

0.3 

0 

0.1 

327.86885 

1 

200 


Table 3: MLE Results from Simulations 


Pulsar 

Type of Run 

Average 9 0 

Percent Error 

Crab 

Stationary 

0.500145 

0.029 

B1821-24 

Stationary 

0.499825 

0.035 

Crab 

Constant V 

0.50012 

0.024 

B1821-24 

Constant V 

0.500745 

0.149 


For the time-varying frequency cases, with constant acceleration, the TOAs arriving at the spacecraft detector 
were broken up into blocks of size T = 1 second. The MLE-DPLL cascade was then applied to track the phase. The 
loop bandwith was taken to be B ^ = 0.05 Hz and the damping ratio was set at £ = 0.707, which is the standard 
underdamped response. The gain coefficients K\ and K 2 for the second-order DPLL were chosen from Table VI in 
[38]. This resulted in I\\ = 0.1199 and K 2 = 0.007658. The loop phase and frequency states were initiated using 
9° pll (0) = 9™ LE ( 0) and f° PLL ( 0) = f s v(t 0 )/c. 

Figure 11 contains the results for the simulated phase tracking using the Crab pulsar with a constant acceleration 
as given in Table 2. The top plots show the phase outputs of the MLE algorithm and the MLE-DPLL algorithm 
along with the truth trajectory. It can be seen that they both track the true phase over the observation window, 
but the error is smaller for the DPLL phase as expected. The bottom two plots show the DPLL frequency state 
versus the true frequency evolution and the error. The DPLL frequency is tracking with some fluxuations, but with 
an error on the order of 10 -4 . Figure 12 contains the results for the simulated phase tracking using pulsar B1821-24 
with a constant acceleration as given in Table 2. The top plots show the phase outputs of the MLE algorithm and 
the MLE-DPLL algorithm along with the truth trajectory. It can be seen that they both track the true phase over 
the observation window, but the error grows linearly as the observation time increases. The error is reduced in the 
DPLL output over the MLE output as expected, but both errors grow without an apparant bound. This might 


17 





be due to the non-zero acceleration and a third order DPLL might need to be considered in future research if the 
phase error is getting out of hand. The bottom two plots show the DPLL frequency state versus the true frequency 
evolution and the error. The DPLL frequency is tracking with some fluxuations, but with an error on the order of 
1CT 3 . 

7.2 Experimental Results 

The same pulsars were used for the experimental runs, using the GXLT, as were used in the computer simulations. 
The effective source and background fluxes are different and the frequency for B1821-24 is slightly different in the 
model used in the GXLT versus the computer simulations. These parameters are all included in Table 4. An amplified 
version of B1821-24 was again used in order to implement the phase tracking algorithm for time- varying frequency. 


Table 4: Pulsar Parameters for Experimental Runs 


Pulsar 

fs (Hz) 

alpha (ph/s 

beta (ph/s) 

B0531+31 (Crab) 

29.9401 

190 

1510 

B1821-24 

327.405595 

0.083 

0.41 

B1821-24 amplified 

327.405595 

1500 

100 


Experimental runs were designed around the same basic three scenarios as with the computer simulations. Table 
5 contains a list of descriptions of all the experimental runs performed. Once again for each of the constant signal 
frequency cases, ten seperate runs were completed and their results were combined with Monte-Carlo methods to 
reduce the chance the performance was an outlier due to non-deterministic noise. The truth trajectories for each 
constant acceleration case were the same as they were for the computer simulations. One extra experiment was 
carried out using the Crab pulsar with a detector moving at a constant acceleration of 5 km/ s 2 with an initial 
velocity of 1 krn/s and an initial phase of 0.5. The goal was to determine if the phase-tracking algorithm could still 
adequatly track the trajectory even with a significantly higher acceleration. This was tested using the crab pulsar 
due to it’s higher flux than B1821-24. 

The results obtained from using the MLE to estimate the initial phase for the stationary and constant velocity 
cases are contained in Table 6. The results all had less than 2.5 % error from the true initial phase of 9q = 0.5. The 
spacecraft TOAs obtained in the constant acceleration case were processed, as in the computer simulations, by the 
MLE-DPLL cascade with a block size of T = 1 second. The gains and parameters for the DPLL were all taken to 
be the same as in the simulations. Figure 13 contains the results for the experimental phase tracking using the Crab 
pulsar with a constant acceleration as given in Table 5. The top plots show the phase outputs of the MLE algorithm 



Figure 10: This is an example of an MLE plot for the crab pulsar with constant velocity v = 1 km/s and initial 
phase 9 0 = 0.5. It is seen that the peak is near 0.5 as desired and the MLE gave an output of 9 0 = 0.50035 in this 
case. 
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Figure 11: This figure shows the results of the simulation using the Crab Pulsar with a detector trajectory along 
the line-of-sight directon with an initial acceleration of 0.1 krn/s 2 , zero inital velocity, and 0.3 initial phase. Plot a 
shows the phase over the 200 second observation. Plot b shows the phase error. Plot c contains the DPLL Doppler 
frequency as it tracks the desired frequency. Plot d shows the error in the DPLL frequency. 


and the MLE-DPLL algorithm along with the truth trajectory. It can be seen that they both track the true phase 
over the observation window, but the error is smaller for the DPLL phase as expected. The phase tracks exactly 
one cycle out of phase, but still tracks the fractional part of the phase accurately, which is the goal. The bottom 
two plots show the DPLL frequency state versus the true frequency evolution and the error. The DPLL frequency 
is tracking with some initial fluxuations, but flattens out to a negligible amount. Figure 14 contains the results for 
the simulated phase tracking using the pulsar B1821-24 with a constant acceleration as given in Table 5. The top 
plots show the phase outputs of the MLE algorithm and the MLE-DPLL algorithm along with the truth trajectory. 
It can be seen that they both track the true phase over the observation window and that the error decreases down 
more for the DPLL phase but takes longer to settle. The bottom two plots show the DPLL frequency state versus 
the true frequency evolution and the error. The DPLL frequency is tracking with some fluxuations, but with an 
error on the order of 0.005 after an initial spike. For the case of a greater constant acceleration and nonzero initial 
velocity as given in Table 5, the results can be seen in Figure 15. The phase is tracking a few cycles out of phase 
this time, but the DPLL phase tracks the fractional part with decreasing error over the observation interval. The 
MLE output is growing as time progresses. It can also be seen that the Doppler frequency is being tracked well by 
the DPLL following an initial fluctuation of over 50 seconds. 


Table 5: Descriptions of the Performed Experimental Runs 


Pulsar 

Type of Run 

th_0 

V0 (km/s) 

a (km/s2) 

Initial f_d (hz) 

Number of Runs 

Observation Time (sec) 

Crab 

Stationary 

0.5 

0 

0 

29.9401 

10 

10 

B1821-24 

Stationary 

0.5 

0 

0 

327.4056 

10 

10 

Crab 

Constant v 

0.5 

1 

0 

29.9402 

10 

10 

B1821-24 

Constant v 

0.5 

1 

0 

327.4067 

10 

10 

Crab 

Constant a 

0.3 

0 

0.1 

29.9401 

1 

200 

B1821-24 

Constant a 

0.3 

0 

0.1 

327.4056 

1 

200 

crab 

Constant a 

0.5 

1 

5 

29.9402 

1 

250 
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Table 6: MLE Results from Experiments 


Pulsar 

Type of Run 

Average thetaJ) 

Error 

Crab 

Stationary 

0.512305 

2.461 

B1821-24 

Stationary 

0.500035 

0.007 

Crab 

Constant V 

0.50816 

1.632 

B1821-24 

Constant V 

0.49387 

1.226 


a 


b 




c xIO ' 3 d 




Figure 12: This figure shows the results of the experimental run using the B1821-24 pulsar with a detector trajectory 
along the line-of-sight directon with an initial acceleration of 0.1 krn/s , zero inital velocity, and 0.3 initial phase. 
Plot a shows the phase over the 200 second observation. Plot b shows the phase error. Plot c contains the DPLL 
Doppler frequency as it tracks the desired frequency. Plot d shows the error in the DPLL frequency. 
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Figure 13: This figure shows the results of the experimental run using the Crab Pulsar with a detector trajectory 
along the line-of-sight directon with an initial acceleration of 0.1 krn/s 2 , zero inital velocity, and 0.3 initial phase. 
Plot a shows the phase over the 200 second observation. Plot b shows the phase error. Plot c contains the DPLL 
Doppler frequency as it tracks the desired frequency. Plot d shows the error in the DPLL frequency. 




Figure 14: This figure shows the results of the experimental run using the pulsar B1821-24 with a detector trajectory 
along the line-of-sight directon with an initial acceleration of 0.1 krn/s, zero inital velocity, and 0.3 initial phase. 
Plot a shows the phase over the 200 second observation. Plot b shows the phase error. Plot c contains the DPLL 
Doppler frequency as it tracks the desired frequency. Plot d shows the error in the DPLL frequency. 
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Figure 15: This figure shows the results of the experimental run using the Crab Pulsar with a detector trajectory 
along the line-of-sight directon with an initial acceleration of 5 km/s 2 , an inital velocity of 1 km/s, and 0.5 initial 
phase. Plot a shows the phase over the 250 second observation. Plot b shows the phase error. Plot c contains the 
DPLL Doppler frequency as it tracks the desired frequency. Plot cl shows the error in the DPLL frequency. 


8 Conclusion 

X-ray navigation provides many opportunities for use in spacecraft applications, either as an autonomous deep 
space navigation system or as a means of reducing the errors as part of a larger filter also receiving DSN telemetry. 
Phase tracking is a method of position determination using the periodic nature of a pulsar’s signal. It has the benefit 
of not relying on time transfer of TOAs to the SSB to compare to a template, and it allows for near continuous 
updates with small blocks of TOAs. In order to be implemented, though, a means of accurately estimating the phase 
and tracking the signal pulses is needed. For a constant frequency observed signal, an MLE is all that is required 
to track the phase. When the dynamical motion of the spacecraft causes the observed signal frequency to vary 
apreciably a cascaded MLE-DPLL algorithm was used where the observation interval was broken up into one second 
blocks. In this paper, these algorithms were tested for a range of scenarios with both the Crab pulsar and B1821-24. 
Simulated results were generated first to confirm the theory of the proposed algorithm and then experiments were 
run using the GXLT equipment at the GSFC. This allowed for actual X-ray photon generation using an MXS and 
detection using a SDD. This added complexity allowed us to test the algorithms with actual photon processing with 
hardware in the loop. In order to test the MLE under constant signal frequency a stationay case and a constant 
velocity case were considered for both pulsars. To test the MLE-DPLL cascade under the varying frequency model, 
a scenario with constant acceleration was considered. 

In Summary, this paper has demonstrated: 

• MLE was able to track the initial phase for either a stationary detector or one moving at constant velocity 
with error under 0.15% in each simulated case, and tracked with error under 2.5% for all of the experimental 
runs . Where all of these values were averaged over ten runs to offset the non-deterministic effects. 

• MLE-DPLL phase tracking algorithm was validated for two experimental pulse sources for constant acceleration 
in one-dimension. 

• MLE-DPLL was shown to accurately track the doppler frequency and the fractional component of the phase. 

• MLE-DPLL tested at an acceleration of 5 km/s 2 with good results when the Crab Pulsar is the source. 
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In the simulations the phase and frequency were tracked for both pulsars. For the Crab pulsar the DPLL phase 
error was on the order of 10~ 3 and doppler frequency error was on the order of 10“ 4 . For B1821-24 the DPLL phase 
error grew to 0.02 cycles error in 200 seconds, which was still better than the MLE alone. The Doppler frequency was 
tracked with frequency error on the order of 10~ 3 . In the experimental runs the phase for the Crab pulsar tracked one 
cycle out of sync, but the phase error settled down to on the order of 10 -2 . The doppler frequency also tracked after 
a bump over the first 50 seconds. For B1821-24 the phase tracked with error on the order of 0.05 and the doppler 
frequency tracked with error on the order of 0.005. An extra scenario with a higher acceleration and a nonzero 
initial velocity was performed experimentally with the Crab. This was to test the algorithm’s durability to higher 
acceleration. The phase tracked with a three cycle offset with the error in the DPLL phase decreasing to 0 over 250 
seconds, but the MLE phase outputs increased over the observation. The doppler phase of the DPLL tracked well 
after a initial transient period of approximately 75 seconds. In all cases the DPLL state estimates were better in the 
long run than just the MLE outputs by themselves. This shows the utility of adding the DPLL in the loop in order to 
decrease the phase error and produce an estimate of the doppler frequency. The results expirementally validated the 
phase-tracking algorithm using lab equipment in place of an actual pulsar and one dimensional trajectories. Future 
work includes looking into methods of speeding up the MLE search algorithm including using fourier transorfsorms 
on the incoming data. Also, DPLLs of higher order will be considered in order to see if they can better deal with the 
effects of nonzero acceleration in some of the test cases. Three dimensional orbits will be another logical experimental 
extension of this work. Overall, the pulse phase tracking algorithm aclieived experimental results matching up with 
the simulated functionality of the algorithm and it would be ideal for use in an XNAV system. 
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